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Abstract 



From extensive molecular dynamics simulations on immiscible two-phase 
flows, we find the relative slipping between the fluids and the solid wall ev- 
erywhere to follow the generalized Navier boundary condition, in which the 
amount of slipping is proportional to the sum of tangential viscous stress 
and the uncompensated Young stress. The latter arises from the deviation 
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of the fluid-fluid interface from its static configuration. We give a continuum 
formulation of the immiscible flow hydrodynamics, comprising the general- 
ized Navier boundary condition, the Navier-Stokes equation, and the Cahn- 
Hilliard interfacial free energy. Our hydrodynamic model yields interfacial and 
velocity proflles matching those from the molecular dynamics simulations at 
the molecular-scale vicinity of the contact line. In particular, the behavior at 
high capillary numbers, leading to the breakup of the fluid-fluid interface, is 
accurately predicted. 
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I. INTRODUCTION 



Immiscible two- phase flow in the vicinity of the contact line (CL), where the fluid- fluid 
interface intersects the solid wall, is a classical problem that falls beyond the framework 
of conventional hydrodynamics [1-12]. In particular, molecular dynamics (MD) studies 
have shown relative slipping between the fluids and the wall, in violation of the no-slip 
boundary condition [6,7]. There have been numerous ad-hoc models [1,8,10-12] to address 
this phenomenon, but none was able to give a quantitative account of the MD slip velocity 
profile in the molecular-scale vicinity of the CL. While away from the moving CL the small 
amount of relative slipping was found to follow the Navier boundary condition (NBC) [13], 
i.e., relative slipping proportional to the tangential viscous stress, in the molecular-scale 
vicinity of the CL the NBC failed totally to account for the near-complete shp. This failure 
casts doubts on the general apphcability of the NBC to immiscible flows and hinders a 
continuum formulation of the hydrodynamics in the CL region. In particular, a (possible) 
breakdown in the hydrodynamic description for the molecular-scale CL region has been 
suggested [7] . In another approach [14] , it was shown the MD results can be reproduced by 
continuum finite element simulations, provided the slip profile extracted from MD is used as 
input. This work demonstrated the feasibility of the hybrid algorithm, but left unresolved the 
problem concerning the boundary condition governing the CL motion. Without a continuum 
hydrodynamic formulation, it becomes difficult or impossible to have realistic simulations 
of micro- or nanofluidics, or of immiscible flows in porous media where the relative wetting 
characteristics, the moving CL dissipation, and behavior over undulating solid surfaces may 
have macroscopic implications. 

Prom MD simulations on immiscible two-phase flows, we report the finding that the gen- 
eralized Navier boundary condition (CNBC) applies for all boundary regions, whereby the 
relative slipping is proportional to the sum of tangential viscous stress and the uncompen- 
sated Young stress. The latter arises from the deviation of the fluid-fluid interface from its 
static configuration [10]. By combining CNBC with the Cahn-Hilliard (CH) hydrodynamic 
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formulation of two-phase flow [11,12], we obtained a consistent, continuum description of 
immiscible flow with material parameters (such as viscosity, interfacial tension, etc) directly 
obtainable from MD simulations. The convective-diffusive dynamics in the vicinity of the 
interface and the moving CL also means the introduction of two phenomenological dynamic 
parameters whose values can be fixed by comparison with one MD flow profile. Once the 
parameter values are determined from MD simulations, our continuum hydrodynamics can 
yield predictions matching those from MD simulations (for different Couette and Poiseuille 
flows). Our findings suggest the no-slip boundary condition to be an approximation to the 
GNBC, accurate for most macroscopic flows but failing in immiscible flows. These results 
open the door to efficient simulations of nano- or microfluidics involving immiscible com- 
ponents, as well as to macroscopic immiscible flow calculations, e.g., in porous media, that 
are physically meaningful at the molecular level [15]. The latter is possible, for example, 
by employing the adaptive method based on the iterative grid redistribution introduced in 
Ref. [15]. This method has demonstrated the capability of resolving, at the same time, both 
the global behavior of a partial differential equation solution with coarse mesh, and a strong 
singularity in a locahzed region with a reflned local mesh of over 10^ ratio to the coarse 
mesh. 



II. MOLECULAR DYNAMICS SIMULATIONS 

The MD simulations were performed for both the static and dynamic conflgurations in 
Couette and Poiseuille flows. The two immiscible fluids were conflned between two parallel 
walls separated along the z direction, with the fluid-solid boundaries defined by 2; = 0, 
(see Fig. 1 for Couette geometry). Interaction between the fluid molecules was modeled 
by a modifled Lennard- Jones (LJ) potential [/// = 4e (cr/r)^^ — 5// (u/r)^ , where r is 
the distance between the molecules, e and a are the energy scale and range of interaction, 
respectively, and 5// = 1 for like molecules and 5// = —1 for molecules of different species. 
Each of the two walls was constructed by two (or more) [001] planes of an fee lattice (see 
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Appendix A), with each wall molecule attached to the lattice site by a harmonic spring. 
The mean-squared displacement of wall molecules was controlled to obey the Lindemann 
criterion. The wall-fluid interaction was also modeled by a LJ potential U^f, with energy 
and range parameters Cwf — 1.16e and awf — l.OAa, and a 5u,f for specifying the wetting 
property of the fluid. Both Uff and U^f were cut off at Tc = 2.5a. The mass of the wall 
molecule was set equal to that of the fluid molecule m, and the average number densities for 
the fluids and wall were set at p = O.Sl/cr^ and — 1.86/(7^, respectively. The temperature 
was controlled at 2.8e/kB, where ks is Boltzmann's constant. Moving the top and bottom 
walls at a constant speed V in the ±x directions, respectively, induced the Couette flow [7] . 
Applying a body force mgext to each fluid molecule in the x direction induced the Poiseuille 
flow [6]. Periodic boundary conditions were imposed on the x and y boundaries of the 
sample. Most of our MD simulations were carried out on samples consisting 6144 atoms for 
each fluid and 2880 atoms for each wall. The sample is 163.5(7 by 6.8(7 along the x and y, 
respectively, and H — 13.6a. Our MD results represent time averages over 20 to 40 miUion 
time steps. For technical details of our MD simulations, we followed those described in Ref. 
[16]. 

Two different cases were considered in our simulations. The symmetric case refers to 
identical wall-fluid interactions for the two fluids (both 6u,f = 1), which leads to a flat static 
interface in the yz plane with a 90° contact angle. The asymmetric case refers to different 
wall-fluid interactions, with S^f — 1 for one and S^f — 0.7 for the other. The resulting 
static interface is a circular arc with a 64° contact angle. We measured six quantities in 
the Couette-flow steady states of = 0.25(e/m)^/^, H — 13.6(7 for the symmetric case and 
V — 0.2(e/m)^/^, H — 13.6a for the asymmetric case: I'j'*^, the slip velocity relative to 
the moving wall; G^, the tangential force per unit area exerted by the wall; the a^x, (^nx 
components of the fluid stress tensor [n denotes the outward surface normal), and v^-, v^. 

We denote the region within 0.85a — zq of the wall the boundary layer (BL). It must 
be thin enough to render sufficient precision for measuring w^'*^, while thick enough to 
fully account for the tangential wall-fluid interaction force, due to the flnite range of the 
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LJ interaction. Thus it is not possible to do MD measurements strictly at the fluid-solid 
boundary, not only because of poor statistics, but also because of this intrinsic limitation. 
The wall force can be identified by separating the force on each fluid molecule into wall- 
fluid and fluid-fluid components. For < z < zq the fluid molecules can detect the atomic 
structure of the wall. When coupled with kinetic coUisions with the wall molecules, there 
arises a nonzero tangential wall force that varies along the z direction and saturates at z zq. 

is the saturated total tangential wall force per unit wall area (Fig. 2). In Appendix A 
we give account of our MD results on both the tangential and normal components of the 
wall force, plus the effect (s) of increasing the wall thickness in our simulations from 2 layers 
of wall molecules to 4 layers and to infinite layers (by using the continuum approximation 
beyond the 4 layers). 

Spatial resolution along the x and z directions was achieved by evenly dividing the 
samphng region into bins, each Ax = 0.425(7 by = 0.85(7 in size, i'^'**' was obtained as the 
time average of fluid molecules' velocities inside the BL, measured with respect to the moving 
wall; was obtained from the time average of the total tangential wall force experienced by 
the fluid molecules in the BL, divided by the bin area in the xy plane; (Txx(nx) was obtained 
from the time averages of the kinetic momentum transfer plus the fluid-fluid interaction 
forces across the constant-a;(2;) bin surfaces, and Vx{z) was measured as the time-averaged 
velocity component (s) within each bin. For the contribution of intermolecular forces to the 
stress, we have directly measured the fluid-fluid interaction forces across bin surfaces instead 
of using the Irving-Kirkwood expression [17], whose validity was noted to be not justified at 
the fluid-fluid or the wall-fluid interface (see the paragraph following equation (5.15) in the 
above reference). In Appendix B we give some details on our MD stress measurements. As 
reference quantities, we also measured cr^^, (t°^ in the static (V = 0) configuration. In 
addition, we measured in both the static and dynamic configurations the average molecular 
densities pi and p2 for the two fluid species in each bin to determine the interface profile. 
The shear viscosity r] = l.QSy^em/cr^ and the interfacial tension 7 = 5.5e/(T^ were also 
determined. 
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We have also measured the interface and velocity profiles for the Poiseuille fiow in the 
asymmetric case, as well as for Couette fiows with different V and H in the symmetric case. 



III. GENERALIZED NAVIER BOUNDARY CONDITION 



In the presence of a fluid-fiuid interface, the static fluid stress tensor (t° contains static 
Young stress (surface tension) as well as those stresses arising from wall-fluid interaction. 
For the consideration of moving CL, we will be concerned with the part of the fluid stress 
tensor which is purely dynamic in origin, i.e., arising purely from the hydrodynamic motion 
of the fluid (and the CL). In the notations below, the over tilde will denote the quantity to be 
the difference between that quantity and its static part. Thus if cr is the total stress, we will 
be concerned only with the hydrodynamic part, denoted hy a — a — cr°. We note that in the 
absence of body forces, force equilibrium in bulk fluid is governed by the relation V • <t = 0. 
As shown in Fig. 2, this relation is altered in the BL, where the fluid-wall interaction means 
the existence of a dynamic, tangential force density g'^ such that V ■ & + g'^x — inside the 
BL. The total tangential force exerted by the wall on the fluid is given by = Jq° dzg"^ per 
unit wall area. In steady state, this wall force is necessarily balanced by the tangential fluid 
force G^. — /q^° dz {dxCTxx + dzdzx)- Here dx,z,n means taking partial derivative with respect 
to X, or surface normal. 

We now present evidences to show that everywhere on the boundaries, relative slipping 
is proportional to G^. (the GNBC, see also equation (3) below): 



where we have used the fact that azx{^) — 0. (More strictly, ^^-^(O ) = because there is no 
fluid below z — Q and hence no momentum transport across z — Qi) Here the z coordinate 



Gi = Pvt 



(1) 



where (5 is the slip coefficient and G| can be written as 




(2) 



7 



is for the lower fluid-solid boundary (same below), with the understanding that the same 
physics holds at the upper boundary, and 9„ — —d^ for the lower boundary 

We have verified the steady state force balance + = on the two boundaries 
(inset to Fig. 3) [18]. This relation reflects the fact that at steady state, the frictional force 
exerted by the solid wall on the moving/slipping fluid is fully accounted for in G^. Thus the 
GNBC (or NBC) can be expressed in either or G!^, but not both. In Fig. 3 we show the 
measured MD data for the symmetric and asymmetric cases in the Couette geometry. The 
symbols represent the values of G^ measured in the BL. The solid lines represent the values 
of Gl calculated from j3v^'^^ by using (3 = Pi = (32 = l.2^/em/a^ for the symmetric case 
and Pi — P2 — 0.532-^/em/(7^ for the asymmetric case away from the CL region 

(straight hne segments in Fig. 3), and P — {Pipi + P2P2)/{pi + P2) in the CL region [19], 
with v^"^^ and pi^2 obtained from MD simulations. It is seen that for the lower boundary 
(upper right panel), the MD data agree well with the predictions of Eq. (1). For the upper 
boundary (lower left panel) the straight line segments also agree well with Eq. (1). However, 
there is some discrepancy in the interfacial region of the upper boundary that seems to arise 
from a "shear thinning" effect of decreasing P at very large tangential stresses [13]. 

The fact that the wall force density is distributed inside a thin BL and vanishes beyond 
the BL necessitates the form of G^ as defined by Eq. (2). However, it is intuitively obvious 
that the fluids would experience almost the identical physical effect (s) from a wall force 
density G'^5{z), concentrated strictly at the fluid-solid boundary with the same total wall 
force per unit area. In the inset to Fig. 2, it is shown that the MD-measured wall force 
density is a sharply peaked function. The sharp boundary limit involves the approximation 
of replacing this peaked function by 5{z). The replacement of a diffuse boundary by a sharp 
boundary can considerably simplify the form of the GNBC, because local force balance 
along X then requires d^d^x + dz^zx — away from the boundary z — Q. Integration of this 
relation from O"'" to zq yields dx Sq° dzdxx{z) + cFzx{zo) — ^zx{^'^) — and as a consequence 
(by comparing with Eq. (2)) G^ = — cr„a;(0"'"). Therefore, a^x changes from 5"2a;(0~) = to 
^zxiS^^) — Gl &X, z — leading ioV-a — GlS{z). Comparing with the diffuse boundary, 
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we see that the form of the equation remains the same, but the BL is now from 0~ to 0"*". 
Thus the GNBC (1) becomes —anxiO) = Pv^''^^ in the sharp boundary hmit. 

The tangential stress dnx can be decomposed into a viscous component and a non- viscous 
component: anxi^) = crnx(-^)+^nx('^)- In Fig. 4 we show that away from the interfacial region 
the tangential viscous stress cr^^i^z) = i]{dnVx + dxVn){z) is the only nonzero component, but 
in the interfacial region a^^ = o'nx — o'nx~'^nx ~ '^nx~'^nx dominant, thereby accounting for 
the failure of NBC to describe the CL motion. Therefore away from the CL region the NBC is 
valid, but in the interfacial region the NBC clearly fails to describe the CL motion. We wish 
to clarify the origin of a^^ and cr^^ as the dynamic and static Young stresses, respectively, 
so that a^^ — a^^ — (7°^ is the uncompensated Young stress. As shown in the inset to Fig. 
4, the integrals (across the interface) of a^^ (= anx — <^nxi calculated by subtracting the 
viscous component r){dnVx + dxVn) from the total tangential stress anx) and a^x are equal to 
7C0S^ci and 7 cos ^5, respectively, at different values of z, i.e., — J^^^dxa^xiz) — j cos 9d{z) 
and — J^^^dxa'l^{z) — 'ycosds{z), where 9diz) and 9siz) are respectively the dynamic and 
static interfacial angles at z [20] . Here Jj^^^ dx denotes the integration across the fluid- fluid 
interface along x. These results clearly show the origin of the extra tangential stress in the 
interfacial region to be the interfacial (uncompensated) Young stress. Thus the GNBC is 
given by 

Pvf^ = -a^xiO) = - [vdnVx] (0) - aUO). (3) 

Here only one component of the viscous stress is nonzero, due to Vn = at the bound- 
ary; and —a^xiO) is denoted the uncompensated Young stress, satisfying — Jmt^nxi^)^^ — 
^{cosO^^'^^ — cos 9f^'^^), with O^^^^-^ being a microscopic dynamic (static) contact angle at the 
fluid-solid boundary. The fact that o^nxi^) ~ ^ away from the CL shows that the GNBC 
implies NBC for single phase flows. 

Due to the diffuse nature of the BL in MD simulations, the contact angle ^^"J)"'^ cannot 
be directly measured. Nevertheless, they arc obtainable through extrapolation by using the 
integrated interfacial curvature within the BL. That is, in the sharp boundary limit the force 
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balance in the fluids is expressed by dxd'xx + dn^nx — 0. Integration in z across the BL gives 
dx r dzaxx(z) - a:,{zo) + <,(0) - a^Zo) + a^O) = 0. (4) 
Integration (of Eq. (4) along x) across the fluid-fluid interface then yields 

A [ r dzaxxiz)] - I dxal^{zo) + / dxa^M + l^d - l^s = 0. (5) 

.JO J J int J int 

where A [/q° dzaxx{z)] is the change of the z-integrated axx across the interface, /C^ and /C^ 
denote the dynamic and static 2;-integrated interfacial curvatures: 

JCd = cos 9d{zo) - cos 9^/''^, 

and 

/C, =cose,(^o) -cos^f^-^. 

Here A[Jq° dzaxx{z)], a'^^{zo), 9d{zo), and 9s{zo) are obtainable from MD simulations, 
/Cs ~ ±2zocos9f^'^f /H for the circular static interfaces, while (T^a.(0) = vi^n'^^xK^) may 
be obtained by extrapolating from the values of tangential viscous stress at z = Zq, 2zo, 
and Szq. Therefore the microscopic dynamic contact angle ^^"'^■^ can be obtained from Eq. 
(5) . In Appendix B 3 we give a more detailed account of the relationship between the MD 
measured stresses and the stress components in the continuum hydrodynamics. The above 
extrapolation is based on this correspondence. 

We have measured the z-integrated axx — ^xx ~ cTxx The dominant behavior is 

a sharp drop across the interface, as shown in Fig. 5 for both the symmetric and asymmetric 
cases. The value of 92^^-^ obtained is 88 ± 0.5° for the symmetric case and 63 ± 0.5° for the 
asymmetric case at the lower boundary, and 64.5 ±0.5° at the upper boundary. These values 
are noted to be very close to 91"^^^. Yet the small difference between the dynamic and static 
(microscopic) contact angles is essential in accounting for the near complete slip in the CL 
region. 

In essence, our results show that in the vicinity of the CL, the tangential viscous stress 
—cr'^x as postulated by the NBC can not give rise to the near-complete CL slip without taking 
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into account the tangential Young stress — cr^^ in combination with the gradient of the (BL- 
integrated) normal stress u^x- For the static configuration, the normal stress gradient is 
balanced by the Young stress, leading to the Young's equation. It is only for a moving CL 
that there is a component of the Young stress which is no longer balanced by the normal 
stress gradient, and this uncompensated Young stress is precisely the additional component 
captured by the GNBC but missed by the NBC. 



IV. CONTINUUM HYDRODYNAMIC FORMULATION 

For Eq. (3) to serve as a boundary condition in hydrodynamic calculations, we need to 
derive the local value of the uncompensated Young stress cr^j;(0) from a continuum formu- 
lation of the immiscible fiow hydrodynamics. Such a formulation is important for studying 
the macroscopic implications of moving CL's under scenarios beyond the capability of MD 
simulations. As a first-order approximation, we formulate a hydrodynamic model based on 
the GNBC and the CH free energy functional [21] that has been successful in the calculations 
of fiuid-fiuid interfacial phenomena: 



^ir(V0)^ + /(0) 



(6) 



where (p — {p2 — pi)/ {p2 + Pi), fi,^') = —\'r4''^ + \u(l>'^i and r, u are parameters which can 



be directly obtained from MD simulations through the interface profile thickness ^ = \J K/r 
[22], the interfacial tension 7 = 2-\/2r^^/3ii, and the two homogeneous equilibrium phases 



given by the condition of df/dcf) — 0, yielding (/)± — ±^r/u (= ±1 in our case). 

To derive the efi^ects of the CH free energy F on immiscible fiow hydrodynamics, let us 
consider a composition field 0(r) . A displacement of the molecules from r to r' = r -|- u(r) 
induces a local change of 0, Scf) — — u • V0, to the first order in u. The associated change in 
F is given by the sum of a body term and a surface term: 

SF = - I dr[g-u]+ I ds [a^,Ui] , (7) 
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where g = jiVcf) is the capillary force density, with ji — 5F/5(f) — —KV'^cf) — r0 + ucf)^ being 
the chemical potential, and 

= -Kdn<l>d,<f>, (8) 

is the tangential Young stress due to the spatial variation of at the fluid-solid boundary 
(i-Ln). Hence the two coupled equations of motion are the Navier-Stokes equation (with the 
addition of the capillary force density) and the convection-diffusion equation for 4>{r): 



Pn 



^^+(v.V)v 



-Vp + V • O-'' + /XV0 + Pmgext, (9) 



^ + v.V0 = MVV, (10) 

together with the incompressibility condition V • v = 0. Here pm is the fluid mass density, p 
is the pressure, cr^ denotes the viscous part of the stress tensor, PmSext is the external body 
force density (for Poiseuille flows), and M is the phenomenological mobility coefficient. 

Four boundary conditions are required to solve Eqs. (9) and (10). Two are given by the 
impermeability condition, i.e., the normal components of the fluid velocity and diffusive flux 
are zero: Vn — and dnP — 0. The form of the other two differential boundary conditions 
may be obtained from the total free energy 

Ft,t[(f)]^F[<l)] + I ds^^f{4>), (11) 

plus our knowledge of GNBC. Here Jwf{4>) is the interfacial free energy per unit area at the 
fluid-sohd boundary. We use ^wf{(l>) — (^7w//2) sin(7r0/2) to denote a smooth interpolation 
between ±A7^//2, with Aj^f = 7m,/(0+) — 7^/(0_) given by — 7 cos 6'|"'''^ (Young's equa- 
tion) . It should be noted that the form of the smooth interpolation has very little effect on 
the final results. Hence we have chosen a simple interpolation function. Similar to Eq. (7), 
the change in Ftot due to the displacement of the molecules from r to r' = r -|- u(r) is given 

by 

SFtot = - j dr[g-n] + j ds [a^iU^^ , (12) 
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where 



Kdn(l) + 



OA (13) 



is the uncompensated Young stress [12] (see below). The continuum (differential) form of 
the GNBC (3) is therefore given by 

fJvf' = -anM = -77 [dnv,] (0) + [L(0)9,0] (0), (14) 

where [L{(j))dx(j)\ (0), with L(0) = Kdncf) + d^wf{(j))/d(f), is the differential expression for 
— cr^^(O) = — a-^j;(0) + a-°3.(0) in equation (3). Here Kdn(t>dx(t> is — c"^a;(0) as seen in Eq. (8), 
and [97tt,/(0)/90]9a;0 = dx'jwfifp) [23] equals to (j°^(0), in accordance with the static force 
balance relation 9^7^/(0) - (7°^(0) = 0. Prom J.^^dx[Kdn(l)dx(f)]{0) = 7C0sef^^ [24] and 
J^^^dxdx^yjf — —^cosOf^'^f, we see that 



/ [L(0)9,0] (0) = 7(cos^r^ - cosC^O, 

J int 



in agreement with [L(0)9j;0] (0) being the uncompensated Young stress. 

Another boundary condition may be inferred from the fact that L(0) = is the Euler- 
Lagrange equation at the fluid-solid boundary for minimizing the total free energy Ftot[4']- 
That is, L(0) = corresponds with the equilibrium (static) condition where 90/9i-|-v- V0 = 
0. The boundary relaxation dynamics of (p is plausibly assumed to be the first-order extension 
of that correspondence for a nonzero L{(j)): 

^ + V • V0 = -r [L(0)] , (15) 

where F is a (positive) phenomenological parameter. 



V. COMPARISON OF MD AND CONTINUUM HYDRODYNAMICS RESULTS 

Motivated by the methods presented in [25,26], a second order scheme is designed to 
solve the CH hydrodynamic model, comprising the dynamic equations and the four bound- 
ary conditions. Details of the numerical algorithm are presented in Appendix C. Besides 
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those parameters which can be directly obtained from MD simulations, M and F are treated 
as fitting parameters, determined by comparison with MD results (values given in the caption 
to Fig. 6). In Figs. 1 and 6 we show that the continuum model can indeed quantitatively re- 
produce the interface and velocity profiles from MD simulations, including the near-complete 
slip {va; in Fig. 6) of the CL, the fine features in the molecular-scale vicinity of the CL, 
and the fast pressure variation in the BL (inset to Fig. 6), with its implied large interfacial 
curvature. We wish to emphasize that for the comparison with the symmetric case, the pa- 
rameters in the continuum model, including those in the GNBC, are directly obtained from 
the MD simulations, whose velocity profiles are then fitted by those from the hydrodynamic 
calculations with optimized M and F values. Thus the comparison with the asymmetric 
(Couette) case, with /32 directly evaluated from MD simulation data, is without adjustable 
parameters. We have also obtained 9'^^^ = 88.1° and 62.8° for the symmetric and asymmet- 
ric (the lower boundary) cases shown in Figs, la and lb, respectively. Both are in excellent 
agreement with their extrapolated values in MD simulations. For the upper boundary in 
the asymmetric case, our calculated 9^^^ = 65.2°, which differs somewhat from the MD 
extrapolation value of 64.5 ± 0.5°. This difference is a reflection of the discrepancy seen in 
Fig. 3. However, it is noteworthy that the difference in the dynamic contact angles does 
not show up in the velocity profiles, which agree well. 

To further verify that the boundary conditions and the parameter values are local prop- 
erties and hence applicable to flows with different macroscopic conditions, we have varied 
the wall speed V. the system size H, and the flow geometry to check that the same set of 
parameters plus the GNBC are valid for reproducing (a) the velocity proflles from a different 
set of Couette-flow simulations in the symmetric conflgurations, shown in Fig. 7, as well 
as (b) the velocity profiles of the Poiseuille fiow simulations in the asymmetric case, shown 
in Fig. 8. The remarkable overall agreement in all cases affirms the validity of GNBC and 
the hydrodynamic model [27], as well as justifies the replacement of the diffuse fluid-sohd 
boundary (force density) by a sharp boundary. 

Another comparison is the dissipation incurred by the moving CL in the Couette flow 
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geometry. Rate of heat generation per unit wall area is given by = Prom 

this we have to subtract a small but constant relative slipping away from the interface, 
^shp ^ 2Vls/ {H + 2/s), where Ig = r]/ (5 is a. shp length for fluid 1 if < and for fluid 2 if 
> 0. The resulting heat generation rate due to the CL is pV'^WsL (for one wall), where 
L is the length of the CL and Wg defines the width of the CL region: 



Thus CL dissipation is equivalent to a segment, ~ H{Ws/ls), of dissipation by single phase 
flow. Figure 9 shows the variation of Ws as a function of capillary number Ca = r]V/'j for 
the symmetric case of Couette flow. Close to Ca ~ 0.1 the value of Wg increases rapidly, 
in good agreement with the MD results, and beyond which the continuum model failed to 
converge. This corresponds to the breakup of the interface observed in MD simulations [28]. 



In summary, we have found for the flrst time the boundary condition that yields near- 
complete slipping of the CL, in good agreement with MD results on the molecular scale. It 
should also be noted, however, that the present continuum formulation can not calculate 
fluctuation effects that are important in MD simulations. Long range interactions, e.g., that 
due to van der Waals interaction, have also been ignored. The latter is potentially important 
in the calculations involving wetting layers. 
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APPENDIX A: WALL-FLUID INTERACTIONS 



We have measured both the tangential and normal components of the wall force exerted 
on the fluids. Both components vary along the z direction and saturate somewhere away 
from the fluid-solid boundary. The tangential component saturates (by 99.8%) at z — zq, 
which is well inside the wall-fluid interaction range {zq — 0.85a, smaller than the cut-off 
distance Tc = 2.5a for the wall- fluid interaction potential Uwf). On the other hand, the 
normal component is 87% of the saturation value at z — zq and 99.8% at ^ = 2zo. The 
different saturation ranges of the tangential and normal components may be understood as 
follows. 

For a fluid molecule close to the solid wall, the interaction with one particular (the 
closest) wall molecule can be much stronger than all the others. As this fluid molecule 
moves laterally but still remaining its close proximity to the wall, it would thus experience a 
strong periodic modulation in its interaction with the wall. This lateral inhomogeneity is an 
important source for the tangential component of the wall force. Away from the fluid-solid 
boundary, each fluid molecule can interact with many wall molecules on a nearly equal basis. 
Thus the modulation amplitude of the wall potential would clearly decrease with increasing 
distance from the wall. Hence the tangential wall force tends to saturate at the relatively 
short range of z c^t Zq. On the contrary, the normal wall force directly arises from the wall- 
fluid interaction, independent of whether the wall potential is "rough" or not. Consequently, 
the normal wall force saturates much slower than the tangential component. 

The MD results presented in this paper were obtained from simulations using solid walls 
constructed by two [001] planes of an fee lattice. We have also carried out MD simulations 
using thicker confining walls. First we changed the number of molecular layers ([001] planes 
of fee lattice) from two to four in constructing each of the two walls. The wall-fluid interac- 
tion potential Uu,f were still cut off at Tc = 2.5(j. It turned out that neither component of the 
wall force shows any noticeable change. The reason is that for the tangential component, the 
two outer planes are too distant to contribute to the roughness of the wall potential, while 
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for the normal component, the fluid molecules closest to the wall are separated from the two 
outer planes by a distance > Tc- Consequently, both the interface and velocity profiles do 
not show any noticeable change. 

Additional wall layers do not contribute to the perceived modulation of the wall potential 
by the fluid molecules. Nevertheless, they can still affect the tangential wall force by modi- 
fying the organization of the fluid molecules near the wall. Such organization is governed by 
the wall-fluid interaction and can be greatly influenced by the normal wall force. To see the 
effects of normal wall force due to additional wall layers, we used four [001] planes of an fee 
lattice plus a half-space continuum in constructing a wall. The flrst four solid layers show the 
atomic structure detectable by the fluid molecules while the half-space continuum models 
the deeper sohd layers. The wall-fluid interaction was modelled as follows. For an in-range 
pair of fluid and wall molecules separated by a distance r < Tc, the interaction potential is 
still Uwf. Here the wall molecule must be from one of the four solid layers. In addition to this 
short-range interaction, the fluid molecules can also experience the long-range interaction 
potential due to (1) the distant wall molecules in the four solid layers and (2) the continuum. 
For (1) we integrated the term in Uyjf over the out-of-range (r > Tc) area of the solid 
layers while for (2) we integrated the same term over the half-space continuum. According 
to this model, only the in-range (r < Tc) part of the solid wall shows atomic structure to a 
fluid molecule while the out-of-range (r > Tc) part is effectively a half-space continuum. 

We found that the effect of the long-range normal wall force (for Syjf > 0) is to attract 
the fluid molecules to the wall. In fact the average number density in the BL can increase by 
3 — 4% once the long-range force is included. As a result, the slip coefficient /3i(2) increases 
by ~ 5 — 15%. This results in small but visible changes in the interface and velocity profiles. 

These tests have convinced us that by using two [001] planes of an fee lattice to model the 
solid wall, we have captured the dominant wall-fiuid interaction. In fact, using two molecular 
layers to model the solid wall has been extensively practiced in the past MD simulations 
[6,7,13,27,29], although in some instances more molecular layers have also been used [30], 
where the accurate modeling of the normal component of the wall-fiuid interaction force is 
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important. 



APPENDIX B: STRESS MEASUREMENTS IN MD SIMULATIONS 

1. Microscopic formula of Irving and Kirkwood 

Irving and Kirkwood [17] have shown that in the hydrodynamic equation of motion 
(momentum transport), the stress tensor (flux of momentum) may be expressed in terms of 
molecular variables as 



<r(r, t) = <TK{r, t) + au{r, t), 



(Bl) 



where ctk is the kinetic contribution to the stress tensor, given by 



rrii 



Y{r,t) 



P 



rrii 



--V(r,i) 



<^(xi-r),/). 



(B2) 



and (Tu is the contribution of intermolecular forces to the stress tensor, given by 



(Tu{r,t) = --^^^(xi-Xj)Fi/(xi-r),/^. 



(B3) 



Here rrii, Pi, and Xj are respectively the mass, momentum, and position of molecule i, 
V(r, t) is the local average velocity, Fj^ is the force on molecule i due to molecule j, / is the 
probability distribution function 

/(xi,---,x^,Pi,---,pjv,i), 
which satisfies the normalization condition 



j dxi • • • dxivdpi • • • dpjv/ = 1, 



and the Liouville equation 



dl 
dt 



-E 



m 9xj ' dpi 



with U being the potential energy of the system, and (■ ■ •, /) means taking the average for 
a probability distribution function /. 
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Although widely employed in the stress measurements in MD simulations, the above 
expression for (Tjj (Eq. (B3)) represents only the leading term in an asymptotic expan- 
sion, accurate when the interaction range is small compared to the range of hydrodynamic 
variation [17]. This can be seen as follows. 

Consider that all the molecules interact via a pair potential ?7pair(-R) such that the in- 
termolecular force = (R,/ R)U^^^^{R) for Xj = Xj + R. Accordingly, Eq. (B3) can be 
rewritten as 

If R R 

cTu{v,t) ^^j dR^Ul^,{R)p('\r,r + R,t), (B4) 
where p^'^^ is the pair density defined by 

It has been shown (see the appendix in Ref. [17]) that according to the definition that dS-cru 
is the force acting across dS, the full expression for ctu is given by 



1 f RR. r 

(Tu(r, t)^ - dK-—U'JR) / dap^^\r - aR, r - aR + R, t) 

2 J R ^ Uo 



(B5) 



It is readily seen that Eq. (B4) may be obtained from Eq. (B5) by keeping only the lowest 
order term in a Taylor's series in a, i.e., 

p(^) (r - aR, r - aR + R, t) ?s p^^) (r, r + R, . 

That means R • VrP*^^^ (r, r + R, must be neghgible compared with p^^^ (r, r + R, i). Here 
R is on the order of the range of intermolecular force. This approximation, however, can 
not be justified at the fiuid-fiuid or the wall-fiuid interface, where R • VrP*-^'' (r, r + R, t) can 
be comparable in magnitude to p^'^\ 

2. Stress measurement in the boundary layer 



In the study of moving CL, it is of great importance to obtain the correct information 
about stress distributions at both the fluid-fluid and the wall-fluid interfaces. Therefore, 
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we have directly measured the x component of fluid-fluid interaction forces acting across 
the x{z) bin surfaces, in order to obtain the xx{zx) component of <Tu. For example, in 
measuring auzx at a given 2;-direction bin surface, we recorded all the pairs of molecules 
interacting across that surface. Here "interacting across" means that the line connecting a 
pair of molecules intersects the bin surface. For those pairs, we then computed (Tuzx at the 
given bin surface using 



where 5s z is the area of ^-direction bin surface, (i, j) indicate all possible pairs of molecules 
interacting across the bin surface, with molecule i being "inside of z^s^" and molecule j 
being "outside of 7,5 Sz' (molecule i is below molecule j), and Fij^ is the x component of the 
force on molecule i due to molecule j. A schematic illustration is shown in Fig. 10. 

For comparison, we have measured the xx and zx components of (Tu using the discrete 
version of Irving-Kirkwood expression (B3): 



where 5v is the volume of sampling bin, i runs over fluid molecules in the sampling bin, j runs 
over fluid molecules in interaction with molecule i, and (• • •) means taking the time average. 
We found that far from the the fluid-fluid and the wall-fluid interfaces, the results based 
on the Irving-Kirkwood expression agree well with those from direct force measurement, 
whereas near the fluid-fluid or the wall-fluid interface, the two results show appreciable 
differences (up to 50%), especially for the zx component at the fluid-fluid interface. 



We want to note the correspondence between the MD-measured stresses and the con- 
tinuum hydrodynamic stress components. This correspondence is essential to obtaining the 





3. Relation of MD-measured stresses 



to the continuum hydrodynamic stress components 
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microscopic contact angle 0^^^ , defined in the continuum hydrodynamic model but not 
directly measurable in MD simulations. 

The GNBC for the difTusc BL is given by 

d r^o 



(B6) 



which involves only MD measurable quantities. To obtain the contact angle 9^^'^^ from MD 
results, we need to interpret the MD-measured quantities in terms of the various continuum 
variables in the hydrodynamic model. In so doing it is essential to note the following. 

(1) (Txx can be decomposed into a molecular component and a hydrodynamic component: 
= + o^xx' ■ Meanwhile, a^^ can be composed into the same molecular component 
and a hydrostatic component: cr^^ = T^x + c^xx- The molecular component Txx exists even 
if there is no hydrodynamic fluid motion or fluid-fluid interfacial curvature. In particular, 
Txx in the BL depends on the wall-fluid interactions. The change of the BL-integrated Txx 
across the fluid- fluid interface equals the change in the wall- fluid interfacial free energy, i.e.. 



d 



Jo 



dzTxx{z) 



-- A7^j = 7i„/(0+) — 7i„/(0_). On the other hand, the hydrody- 
namic component in (Txx results from the hydrodynamic fluid motion and fluid-fluid 
interfacial curvature. In the static (l^ = or g^xt = 0) configuration, becomes the 
hydrostatic component in a^^. 

(2) (Tzx{zo) can be decomposed into a viscous component plus a Young component: 
(yzx{zo) = crl^izo) + cr^^izo) with = rjid^Vx + dxV^) and /j^t dxal^{zo) = 'y cos9d{zo). 

(3) cr'^xi^o) is the static Young stress: i.e., Jl^^. dxa'^j.{zo) — ^cos9s{zo). 

With the help of the above relations, integration of Eq. (B6) across the fluid-fluid 
interface yields 



/ dxpvfP = A [ dza^f] + [ dxal^izo) + 7 cos ^^(^o) - ^\ T dza. 



HS 

XX 



- 7C0S (2^0), 
(B7) 



where A 



Jq° dza^J^^^^^ is the change of the z-integrated ct^^^^^'> across the interface. 
According to the Laplace's equation, the change of the hydrostatic 2;-integrated normal 
stress is directly related to the static 2;-integrated curvature /C^: 
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-A / c^zfjf/ = 7^. = 7 cos 9,{zo) - cos Of^'^] . (B8) 

Note that /Cs vanishes in the symmetric case. Substituting Eq. (B8) into Eq. (B7) then 
yields 

/ dx(5vf'^ = A r dza^f + / dxa^^izo) + 7 cos ^^(^^o) - 7 cos ^f''-^. (B9) 

J int JQ J int 

If interpreted in the continuum hydrodynamic formulation with a sharp fluid-solid boundary, 
the last term in the right-hand side of Eq. (B9), — 7C0S^|"''-^, is the net wall force along 
X arising from the wall-fluid interfacial free energy jump across the fluid-fluid interface, in 
accordance with the Young's equation cos6'^™'-^ = A^y^f. On the other hand, the sum 
of the first three terms on the right-hand side of Eq. (B9) is the net fluid force along x 
exerted on the three fluid sides of a BL fluid element in the interfacial region, due to the 
hydrodynamic motion of the fluids. 

To obtain an extrapolated value for the contact angle O^^^-^ from Eq. (B9), we turn to 
the Stokes equation in the BL: 

-d,p + d^al, + d,al, + lidA = 0, (BIO) 

obtained from the x-component of Eq. (9) by dropping the inertial and external forces. 

Integration in z of Eq. (BIO) across the BL, together with the integration along x across 
the fluid-fluid interface yields 

A [ r dz {-p + alj] + [ dxal^{zo) + 7cos^<i(^o) - / dxal^{Q) - -fcose^^ = 0. (Bll) 

.Jo J J int J int 

Here we have made use of two relations: (1) iidx(f) — jKd{x — Xi^t) in the sharp interface limit 
[31], with K being the interfacial curvature and Xint the location of the interface along x. (2) 
Jq° dzK is the dynamic 2;-integrated curvature JCd = cos6d{zo) — cos^^^"^-^. The local force 
balance along x is expressed by Eq. (BIO). Accordingly, the force balance along x for the 
BL fluids in the integration region is expressed by Eq. (Bll), where A [Jq° dz {—p + crl^)] is 
the net force on the left and right (constant-x) surfaces, Jint d'Xa^^{zo) + 'ycos9d{zo) is the 
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tangential force on the z — zq surface, and — Jj^^. dxal^{0) — 7 cos is the tangential force 
on the z = surface. Substituting Eq. (Bll) into Eq. (B9) and identifying the normal 
stress —p + cr^j. with a^J^, we obtain 

/ dx(3v^^ = I dxal^{Q) + 7 cos ^^"^ - 7 cos Of^'^ , (B12) 

J mi Jint 

which is identical to the integration of the continuum GNBC (Eq. (14)) along x across the 
fluid-fluid interface. 

In summary, to obtain Eq. (B12) from Eq. (B7), we have used both d^a^f + S^u^^ — 
and d^cr^J^ + dzCr^x = 0, whose integrated expressions are given by Eq. (B8) and Eq. (Bll), 
respectively. We note that dx{cr^J^ — o^^) + dz{(Jzx — cr^zx) = is equivalent to the relation 
dx^xx + dnd'nx — (integrated expressions given by Eqs. (4) and (5)), which has been used 
to obtain 6'^'^'-^ through extrapolation in Sec. III. 



APPENDIX C: NUMERICAL ALGORITHM 



We present our numerical algorithm for solving the continuum hydrodynamic model, 
comprising the dynamic equations (9) and (10) and the four boundary conditions Vn — 
0, dn/jL — 0, plus Eqs. (14) and (15). We pay special attention to the apphcation of 
boundary conditions, and restrict our analysis to the Couette flow because the generalization 
to Poiseuille flow is straightforward. 



1. Dimensionless hydrodynamic equations 

To obtain a set of dimensionless equations suitable for numerical computations, we scale 



4> by \4>±\ = \frju-i length by ^ = ^jKjr^ velocity by the wall speed V ^ time by ^/V^, and 
pressure/stress by r]V/$,. In dimensionless forms, the convection-diffusion equation reads 

^ + V V0 = Ca^\-V^ct> - (/) + 0^), (CI) 

the Navier-Stokes equation reads 
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n 



'dt 



+ (v • V) V 



-Vp + + i3(-VV - + 0'^) V0, 



(C2) 



the relaxation of (p at the fluid-sohd boundary is governed by 



9n0- -g-COSe°S^(0) 



(C3) 



and the GNBC becomes 



■\/2 

dn(t> - — cos6l°s^(0) 



dx(t>- dnVx 



(C4) 



Here = (1 - (fyPi/^ + (1 + 0)/92/2 and 5^(0) = (7r/2) cos(7r(/)/2). Five dimensionless 
parameters appear in the above equations. They are (1) Cd — Mr/V^, which is the ratio 
of a diffusion length Mr/V to ^, (2) TZ = pV^/r], (3) B = r^^/uriV = 3'y/2^/2riV, which is 
inversely proportional to the capillary number Ca — rjV/^, (4) = KT/V, which is the 
ratio of KT (of velocity dimension) to V, and (5) £s(0) = v/Pi'f)^: which is the ratio of the 
slip length ls{(j)) = v/Pi'P) to 



2. Finite-difference scheme 

For immiscible Couette flows, there are four variables 0, v^, v^, and p to be solved 
in a two-dimensional (2D) system (in the xz plane). We want to solve the convection- 
diffusion equation and the Navier-Stokes equation in a 2D system of length (along x) 
and height (along z). Here must be large enough to allow the single phase flows (far 
from the fluid-fluid interface) to approach uniform shear flows. A flnitc-diffcrcnce scheme 
is employed as follows. (1) and A^^ equally spaced levels are introduced in the x and 
z directions, respectively. Grid size is given by — Lx/{Nx — 1) and = Lz/{Nz — 1) 
along X and z, respectively. (2) Each variable {q) is deflned at x sites distributed 
from X — —Lx/2 to and from z — —Lz/2 to represented by the array g^j, with 

i = 1, Nx and j = 1, N^. Here g^j = q{xi, Zj), with = (i - l)Lx/{Nx - 1) - Lx/2 
and Zj = {j — l)Lz/{Nz — 1) — (3) In applying the various boundary conditions, 

"ghost" sites outside the system, i.e., i — 0, i — + 1, j — 0, or j — + 1, may appear 
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in the discretization scheme. The values of the variables at the ghost sites are determined 
separately from the various boundary conditions, detailed below. (4) First and second 
spatial derivatives along ( {=x or z) are represented by d(^q{(k) = [?(Cfe+i) ~ ?(Cfe-i)]/2A^ 
and d^qia) = [q(Ck+i) + q(Ck-i) - 2g(a)]/A^. 



3. Convection-diffusion equation 



With the chemical potential given by 



1^' 



A2 



+ 



A? 



(C5) 



the discretized convection-diffusion equation is 



^0,,,- + [v • V4>\^ = 



A2 

X 



A2 

z 



(C6) 



with 



[v • V0] . . = v^ij 



2A, ■ '''' 2A, 



(C7) 



The boundary conditions ai x = ±Lx/2 can be easily applied using = ±1 and 



vU 



(C8) 



Lz + 2L5 Lz 

for single-phase uniform shear flows. Here we focus on the boundary conditions aX. z = 
dn/jL — and Eq. (C3). We spell out the numerics for the lower boundary j — 1, 
with the understanding that the same can be applied to the upper boundary. 

To solve the discretized convection-diffusion equation (C6) at the lower boundary j = 1, 
we need the values of /lij at j = 1 and j — 0. We also need the values of /lij at j = 1 
to solve the same equation at j — 2. According to Eq. (C5), /lij at j = 1 and j — can 
not be directly evaluated from (pij with i = I, ■■■,Nx and j = 1, A^^. But they can still 
be determined from the boundary conditions at 2; = —Lz/2. /lij at j — is obtained from 
dn/ji = at j = 1 as 
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fJ'i,j-l=0 — fJ'i,j+l=2- 



(C9) 



To obtain /lij at j — 1, we need to determine (pij at j — 0. This can be done by requiring 
that Eqs. (CI) and (C3) yield the same d(j)/dt at z = —Lz/2. The discretized convection- 
diffusion equation is given by Eq. (C6) while the discretized relaxation equation for at 
the boundary j — lis given by 

d 



m 



+ [v • V0] 



-V. 



2A 



4>i,j+i _ \/2 
3 



(CIO) 



Equating the right-hand side of Eq. (C6) at j — 1 (with //j^o fixed by Eq. (C9) and other 
//'s given by Eq. (C5)) with that of Eq. (CIO) leads to a tridiagonal system of linear 
equations for 4>ij (0jj coupled with (pi-ij and 0i+ij) at j = 0. Solving this tridiagonal 
system determines (pij at j — 0, from which we obtain //j j at j = 1 by using Eq. (C5). 



4. Navier-Stokes equation 

We now turn to the Navier-Stokes equation (C2) with the incompressibility condition 
V • V = 0. The difficulty in solving the Navier-Stokes equation is the lack of a time evolution 
equation for the pressure p. In the following, we will introduce a numerical method based 
on the pressure Poisson equation [25]. 



a. Pressure Poisson equation 

Taking the divergence of momentum equation (C2) and applying the incompressibility 
condition, we obtain the pressure Poisson equation 

V^p = -7^V • [(v • V) v] + i3V • [(-VV - + 0^)V0]. (Cll) 

Dotting the momentum equation (C2) with the surface normal at the fluid-solid boundary 
and using Vn — 0, we obtain for Eq. (Cll) the boundary condition 

dnP = + B(-VV - + 0')5n</', (C12) 
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at z — ±Lz/2. In addition, we use V^p = and dxP — for the values of V^p and dnP at 
the boundaries x — ±Lx/2. This reflects the single phase flow given by Eq. (C8). 

From the momentum equation (C2) and the pressure Poisson equation (Cll), we derive 
a diffusion equation 

K^ = V=(V.v), 

for V • V. With V • v = given at time t — 0, and in order to ensure that v remains 
divergence-free at i > 0, we must impose the additional boundary condition V • v = at 
all times t > 0. We will show that this boundary condition is needed in solving for p in a 
finite-difference scheme. 

In order to solve the pressure Poisson equation, we need to evaluate [V^p]i,j for i — 
1, and j = 1, N^, [da:p]ij for i = l,Nx and j = 1, N^, and [dzp]i,j for i = l,...,Nx 
and j — l,Nz. For V^p, we have 

for i = l,Nx and j = 1, A^^; 

Vxi+l,j Vxi—lJ f^ij-i-i '^zi,j—i '^zi+lj '^zi—l,j '^xi,j+l "^xij—l 



2Ax 2A, 2Ax 2A 



z 



,»„ / - '^<l>i,j + <l>i-l,j , <l>i,j+l - '^(t>i,j + (t>i,j-l \ 

I » [ ~ 4>i+l,j - 4>i-l,j , ~ ~ 4>i,j-l 

\ 2Ax 2Ax 2Az 2Az 

for i = 2, A^^ — 1 and j — 2, A^^ — 1; and 



[V^p]i,j = 27^- 



2A, 2A, 



^ ^^'^ V A2 + A? 



z 



2Ax 2Ax 

for i — 2, A^j. — 1 and j — 1, A^;^ (where = and ^^/x = 0). We see that (f) and at 
ghost sites of j — Q^N^ appear in the last expression. The ghost 0's have already been 
determined in solving the convection-diffusion equation, while the ghost v^'s are determined 
through the additional boundary condition V • v = 0: 

'^xi+l,j ~ '^xi-l,j . '^zi,j+l ~ '^zi,j—l p. 

2A^ 2A^ 
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for i — 2, Nx — 1, and j — 1, N^. For we have 
for i — l,Nx and j — 1, Nz] 



for i — l,Nx and j = 1, N^; and 



for i = 2, A^a: — 1 and j — 1, (where ^2 = 0). The last expression involves the ghost 0's 
and v^'s at j — 0,Nz + 1. Given the above values of [V^pjjj and [dnp]i,j, we apply a 2D Fast 
Fourier Transformation to solve Pij{0) (up to a constant) for i = 1, A^^. and j = 1, ...,Nz. 

b. Slip boundary condition 
The discretized Navier-Stokes equation is given by 

~dr ~ 2A; 2A; ^ 2Ax 

^7^^''^■ 2A, 
for i — 2, A^a; — 1 and j — 1, AT^, and 

dVzi,j _ "f^zi+l J ~ "^zi-lj VziJ+1 — Vzi,j-1 1 Pi,j+l — Pi,]-! 

"dT ~ 2K'x Wz n Wz 

1 / 'Vzi+l,j — ^Vzij + Vzi-l,j Vzi,j+1 — 2Vzi,j + Vzi,j-l\ . 

I ^ Al J ^^^^^ 



7^^''^ 2A, 

for i = 2, A^a; — 1 and J — 2, Nz — 1, together with the boundary conditions that Vzij = 
at j = 1, A^2 and v given by Eq. (C8) at i = 1, A^-^. Equation (C13) at j = 1, A^^: involves (f) 
and Vx at ghost sites of j = 0, A^^ + 1. The ghost 0's come from at j = 1, Nz, and have 
already been determined. The ghost Vx^s are determined from the discretized GNBC 
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2A, 



2A, 



(C15) 



at the lower boundary j — 1 with v^^lj — v^ij — V, and 



h,j+i - V2 



2A, 



0, 



i+lj '^xi,j+l ^xi,j—l 



2A, 



2A, 



(C16) 



at the upper boundary j — with v^l^ — Vxij + V. 

In summary, to solve the dynamic equations (9) and (10), we need to use cf) — ±1 and 
Eq. (C8) at a; = ±^3-/2, with f„ = 0, d^H = 0, plus Eqs. (14) and (15) at 2; = ±LJ2. In 
particular, in applying the boundary conditions at 2; = ±Lz/2, values of 0, v^, and at 
ghost sites have to be introduced and solved for. 



5. Time integration 

We outline the scheme for time discretization and integration. For simplicity we only de- 
scribe the forward Euler time stepping. In the following a superscript n denotes consecutive 
time instants and At is the time interval. 

Time Stepping: Given {^i^j} and {v^^ j at all the sites (i = 1, A^^ and j — 1, A^^) in 
the system: 

Step 1: Determine {^i^j}; and {v^jj at the ghost sites from the various boundary 

conditions, as described in Sees. C3, C4a, and C4b. 

Step 2: Solve at all the interior sites (i — 1, and j — 1, N^) from Eq. (Cll) 

with appropriate boundary conditions for dnP, as described in Sec. C4a. 
Step 3: Compute j^i^/^} and jv^^^j at all the interior sites (except those fixed by the 
boundary conditions at all times) using 
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and 



n+l _ n 

n — — — = -n (v" • V) v" - Vp" + vV" + ^/x"V(/.", 



according to Eqs. (C6), (C13), and (C14) in discretized time. Here the ghost |a*"j }, {^ilj }) 
and |v"j | determined in Step 1 and {p^^j} solved in Step 2 are needed. 
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FIGURES 
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FIG. 1. Segments of the MD simulation sample for the immiscible Couette flows. The colored 
dots indicate the instantaneous molecular positions of the two fluids projected onto the xz plane. 
The black/gray circles denote the wall molecules. The upper panel illustrates the symmetric case; 
the lower panel illustrates the asymmetric case. The red circles and the blue squares represent the 
time-averaged interface profiles, defined by pi = P2 {4> = 0), for the two cases. The black solid lines 
are the interface profiles calculated from the continuum hydrodynamic model with the GNBC. 
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FIG. 2. By subdividing the boundary layer into thin sections, we plot the accumulated wall 
force per unit wall area as a function of distance z away from the boundary. Here {£) is defined 
by G'^{z) = Jq dz'g^{z') where g'^ is the density of tangential wall force. For different x positions, 
the absolute value of the saturating total wall force is different. However, when normalized by the 
corresponding saturated total wall force per unit area at each x, all points fall on a universal curve, 
nearly independent of x. It is seen that at z = the wall force has reached its saturation value. 
Inset: Tangential wall force density plotted as a function of distance z away from the boundary. 
The solid lines are averaged g'^ in thin sections at different x, normalized by the corresponding 
saturated total wall force per unit area. The dashed line is a smooth Gaussian fit. In the sharp 
boundary limit this peaked wall force density is approximated by G'^5{z). 
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FIG. 3. (3iV/Gl plotted as a function of V/v^^^^. Symbols are MD data measured in the BL 
at different x locations, where the red circles denote the symmetric case and the blue squares the 
asymmetric case. The solid lines were calculated from Eq. (1) with values of /3i^2 and the expression 
of P given in the text. The statistical errors of the MD data are about the size of the symbols. 
The upper-right data segment corresponds to the lower boundary, whereas the lower-left segment 
corresponds to the upper boundary. The slopes of the two dashed lines are given by /3f2- Inset: 
plotted as a function of , measured in the two BL's at different values of x. The symbols 
have the same correspondence as in the main figure. The data are seen to lie on a straight line 
with a slope of —1, indicating -|- = 0. 
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FIG. 4. Two components of the dynamic tangential stress at z = zq, plotted as a function of 
X. The dashed lines denote dj^; solid lines represent the viscous component. Here red indicates 
the symmetric case and blue the asymmetric case. In the CL region the non-viscous component 
is one order of magnitude larger than the viscous component. The difference between the two 
components, however, diminishes towards the boundary, z = 0, due to the large interfacial pressure 
drop (implying a large curvature) in the BL, thereby pulling 6ij_ closer to Og. Inset: 'Sd,s plotted 
as a function of 7cos0d^s at different values of z. Here 11^ = —/ dx{anx — f^nx)' = — J dxa^^, 
and O^^s was measured from the time-averaged interfacial profiles (Fig. 1). The red circles denote 
the symmetric case, the blue squares the asymmetric case, the solid blue squares the asymmetric 
static case, and the single solid red circle at the origin denotes the symmetric static case. The data 
are seen to follow a straight (dashed) line with slope 1, indicating T,d,s = "J cos 9d,s- 
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FIG. 5. S = axx{z)dz = j^^ [crxx{z) — ^"^^(z)] dz plotted as a function of x. Here red circles 
denote the symmetric case and blue squares the asymmetric case. For clarity, vertically 
displaced such that a^x = far from the interface in the symmetric case, and for the asymmetric 



at the center of the interface. 
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FIG. 6. Comparisons between the MD (symbols) and the continuum hydrodynamics (solid 
lines) results for the Couette flow, the latter calculated with the GNBC and values of 
M = 0.023a'^ /^/me and F = 0.6Qa/^/me. (a) The Vx profiles for the symmetric case 
{V = 0.25(e/m)^/^ and H = 13.60") at different z planes. The profiles are symmetric about the 
center plane, hence only the lower half is shown for z = 0.425fT (black circles), 2.125(7 (red squares), 
3.825(7 (green diamonds), and 5.525fT (blue triangles), (b) The Vx profiles for the asymmetric case 
{V = 0.2(e/m)V2 and H = 13.6cr) at z = 0.425(7 (black circles), 2.975(7 (red squares), 5.525(7 
(green diamonds), 8.075o" (blue triangles), 10.625cr (yellow triangles), 13.175cj (maroon triangles). 
For the boundary layers, Vx = means complete slip. Inset: Pressure variation in the BL for the 
symmetric case. The solid line represents the BL-averaged hydrodynamic pressure Zq^ Jq° p{z)dz 
from the continuum model, and red circles denote Zq^ /q^° axx{z)dz measured in MD simulations 
(see Fig. 5). Note the fast variation across the interface. The interfacial pressure drop in the BL 
is a factor 5 — 10 larger than that in the middle of the sample, implying large curvature. 
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FIG. 7. Comparisons between the MD (symbols) and the continuum hydrodynamics (solid 
lines) results for two symmetric cases in the Couette flow geometry. Compared with Fig. 6a, 
V and H have been varied, respectively, but the continuum results are calculated with the same 
set of parameters and the GNBC. The profiles are symmetric about the center plane, hence only 
the lower half is shown, (a) The Vx profiles for V = 0.25(e/m)^/^ and H = 10.2(T, shown at 
z = 0.425(7 (black circles), 2.125(7 (red squares), and 3.825(7 (green diamonds), (b) The Vx profiles 
for V = 0.275(e/m)"'^/^ and H = 13.6(7, shown at z = 0.425(7 (black circles), 2.125(7 (red squares), 
3.825(7 (green diamonds), and 5.525(7 (blue triangles). 
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FIG. 8. Comparisons between the MD (symbols) and the continuum hydrodynamics (solid 
lines) results for an asymmetric case in the Poiseuille flow geometry. Compared with Fig. 6b, 
the type of flow has been changed, but the continuum results are calculated with the same set 
of parameters, (a) A segment of the instantaneous configuration in the MD simulation. The two 
walls, separated hy H = 13. 60", move at a constant speed V = 0.51(e/?7i)^/^ in the —x direction in 
order to maintain a time-independent steady-state interface, with mgext = O.OSe/ci applied in the 
X direction. The symbols have the same correspondence as those in Fig. lb. The black solid line 
is the interface profiles calculated from the continuum hydrodynamic model. The colored dashed 
lines indicate the z coordinates of the Vx profiles shown in (b). (b) The Vx profiles at z = 0.425cr 
(black circles), 2.125(T (red squares), 3.825cr (green diamonds), and 5.525(7 (blue triangles). The 
profiles are symmetric about the center plane, hence only the lower half is shown. 
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FIG. 9. Width for the moving CL region, Wg, plotted as a function of the capillary number 
Co = 'r]V/'j for the symmetric case by varying V and keeping H = 13. 60". We note that for most 
of the MD data measured in the symmetric case, Ca ~ 0.088. Solid line was calculated from the 
immiscible hydrodynamic model employing the GNBC; red circles denote the MD results. 
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FIG. 10. Schematic illustration of measuring the zx component of cru- The horizontal solid 
lines (separated by short vertical lines) represent bin surfaces with surface normals along the z 
direction. Circles denote fluid molecules. The dashed lines connect pairs of interacting molecules. 
Here the bin surfaces and the molecules are projected onto the xz plane. Molecules that appear 
to be close to each other may not be in the interaction range if their distance along y is too large. 
A pair of interacting molecules may act across more than one bin surface. Here the (1,3) pair acts 
across the surfaces A and C while the (1,5) pair acts across the surfaces B and D. At each bin 
surface the stress measurement must run over all the pairs that act across that surface. For surface 
D, there are three pairs of interacting molecules (1,5), (2,4), and (2,5) that contribute to the zx 
component of au- 
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